
library(gnorm)
library(EnvStats)
library(PLCE)
rm(list=ls())
options(device="quartz")

make.rank <- PLCE:::make.rank
pos_measure <- PLCE:::pos_measure
set.seed(1)
n.use<-20
nr.use<-50000
m1<-matrix(rgnorm(n.use*nr.use,beta=1),nrow=nr.use)
m2<-matrix(rgnorm(n.use*nr.use,beta=2),nrow=nr.use)
m3<-(matrix(rgnorm(n.use*nr.use,beta=8),nrow=nr.use))
#m3<-(matrix(-2*runif(1000*500),nrow=nr.use))*.2
#m4<-(matrix(revd(n.use*nr.use),nrow=nr.use))
m4<-(matrix(rchisq(n.use*nr.use,df=5),nrow=nr.use))
m5<-rbind(m3,m3,m3,m3,m1)


m1<-m1/sd(as.vector(m1))
m2<-m2/sd(as.vector(m2))
m3<-m3/sd(as.vector(m3))
m4<-(m4-mean(m4))/(sd(as.vector(m4)))
m5<-(m5-mean(m5))/(sd(as.vector(m5)))


# Supressing warnings b/c function is internal to PLCE and expecting
# a PLCE object

p1<-suppressWarnings(pos_measure(m1,trim=0.0125))
p2<-suppressWarnings(pos_measure(m2,trim=0.0125))
p3<-suppressWarnings(pos_measure(m3,trim=0.0125))
p4<-suppressWarnings(pos_measure(m4,trim=0.0125))
p5<-suppressWarnings(pos_measure(m5,trim=0.0125))

gg_color <- function(n,alpha0=1) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100,alpha=alpha0)[1:n]
}

pdf("SuppMaterialsFig1_Ratkovic.pdf",h=4*1.1,w=7*1.5)

par(mfrow=c(1,2))

cols<-gg_color(5)

par(mar = c(3, 3,1,0.2), # Dist' from plot to side of page
    mgp = c(2, 0.4, 0), # Dist' plot to label
    las = 1, # Rotate y-axis text
    tck = -.01, # Reduce tick length
    xaxs = "i", yaxs = "i", # Remove plot padding
    oma =c(.2,.2,0,0))

plot(density(as.vector(m1)),col="red",xlim=c(-3,3),ylim=c(0,.8),xlab="",ylab="",main="",type="n",xaxt="n",yaxt="n")


rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "gray90")
abline(h=(1:100)/10,col="white")
abline(v=(-90:90)/1,col="white")
axis(1)
axis(2)
mtext(side=2,line=1.7,font=2,text="Density",las=0)
mtext(side=1,line=1.7,font=2,text="Value")


lines(density(as.vector(m1)),col=cols[1],lwd=2)
lines(density(as.vector(m2)),col=cols[2],lwd=2)
lines(density(as.vector(m3)),col=cols[3],lwd=2)
lines(density(as.vector(m4)),col=cols[4],lwd=2)
lines(density(as.vector(m5)),col=cols[5],lwd=2)

legend.txt<-c("Thin Tailed","Normal","Fat Tailed","Skewed","Mixture")
legend("topleft",legend=legend.txt,col=gg_color(5),lty=1,xpd=T,lwd=2,bg="gray90",bty="n",box.col="gray90")



plot(make.rank(p1$diff),p1$diff,type="n",col="red",lwd=2,ylim=range(c(p1$diff,p2$diff,p3$diff,p4$diff)),
      xaxt="n",yaxt="n",xlab="",ylab="",xlim=c(-5,105))
#abline(c(0,1),col="gray50")
abline(h=0,lwd=2)
rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "gray90")
abline(h=(-100:100),col="white")
abline(v=(0:100)*10,col="white")

mtext(side=1,line=1.7,font=2,text="Percentile",las=0)
mtext(side=2,line=1.7,font=2,text="Excess Kurtosis",las=0)
axis(1)
axis(2)
lines(make.rank(p1$diff),p1$diff,col=cols[1],lwd=2)

lines(make.rank(p2$diff),p2$diff,col=cols[2],lwd=2)
lines(make.rank(p3$diff),p3$diff,col=cols[3],lwd=2)
lines(make.rank(p4$diff),p4$diff,col=cols[4],lwd=2)
lines(make.rank(p5$diff),p5$diff,col=cols[5],lwd=2)

legend.txt<-c("Thin Tailed","Normal","Fat Tailed","Skewed","Mixture")
legend("topleft",legend=legend.txt,col=gg_color(5),lty=1,xpd=T,lwd=2,bg="gray90",bty="n",box.col="gray90")

dev.off()


